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Abstract. With the recent detection of cosmic shear, the 
most challenging effect of weak gravitational lensing has 
been observed. The main difficulties for this detection were 
the need for a large amount of high quality data and 
the control of systematics during the gravitational shear 
measurement process, in particular those coming from the 
Point Spread Function anisotropy. 

In this paper we perform detailed simulations with the 
state-of-the-art algorithm developed by Kaiser, Squires 
and Broadhurst (KSB) to measure gravitational shear. We 
show that for realistic PSF profiles the KSB algorithm can 
recover any shear amplitude in the range 0.012 < \^\ < 
0.32 with a relative, systematic error of 10 — 15%. We 
give quantitative limits on the PSF correction method as 
a function of shear strength, object size, signal-to-noise 
and PSF anisotropy amplitude, and we provide an auto- 
matic procedure to get a reliable object catalog for shear 
measurements out of the raw images. 
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1. Introduction 

Weak gravitational lensing has become one of the most 
important cosmological tools in the last decade. The mass 
distribution in the universe can in principle be mapped 
from the measurement of the coherent distortion caused 
by inhomogeneously distributed foreground mass on the 
orientation of the faint background galaxies. Lensing by 
clusters of galaxies, strong enough to be routinely de- 
tected now, already provides important constraints on 
the nature of our universe [e.g. it probed massive high- 
rcdshift clusters (Luppino & Kaiser 1997; LK97 hence- 



Send offprint requests to: erben@mpa-garching.mpg.de 



forth; Clowe et al. 1998) and revealed the possible ex- 
istence of dark clumps (Erben et al. 2000, Umetsu et 
al. 2000)]. Since the first publications of "cluster lens- 
ing" (see e.g., Tyson, Valdes & Wenk 1990 and Kaiser & 
Squires 1993) much theoretical and observational progress 
on weak lensing has been made, and interest progressively 
turned to very weak shear measurements. Galaxy-galaxy 
lensing is an example of such a study which provides con- 
straints on average dark matter halo properties of galax- 
ies, and cosmic shear is, ultimately, the direct measure- 
ment of statistical properties of the large-scale matter 
distribution in our universe (see, e.g., Schneider & Rix 
1997, Schneider et al. 1998b). First galaxy-galaxy lensing 
experiments investigating dark halos of field and cluster 
galaxies have been successfully made (Brainerd, Bland- 
ford & Small 1995, Fischer et al. 1999, Natarayan et al. 
1998). Very recently, the cosmic shear was also detected 
(Schneider et al. 1998a, van Waerbeke et al. 2000, Bacon 
et al. 2000a, Wittmann et al. 2000, Kaiser et al. 2000). 
Because of the small amplitude of the lens effect, a quanti- 
tative analysis of the signal requires high-quality data and 
very precise data analysis for the detection of the signal. 
With the advent of new high-quality wide field CCD cam- 
eras substantial amounts of useable data for these stud- 
ies will soon become available [e.g., the Descart project 
(http://terapix.iap.fr/Descart/Descart_english.html)]. 

On the analysis side, the key issue is the determination 
of the gravitational shear 7 measured from galaxy shapes. 
Using real data, we have to deal with observational and 
data reduction defects. On CCD images, the information 
from very faint and small objects whose shape we want 
to measure is often contained in only a few image pix- 
els and a clear relation of pixels to individual objects is 
hard to achieve. Moreover, the measurements are rendered 
more difficult by pixel noise and PSF effects that can 
mimic a possible shear signal. Kaiser, Squires & Broad- 
hurst (1995, KSB henceforth) have proposed a method for 
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shear determination which is optimised for the analysis of 

faint and small objects by utilising weighted brightness 
moments of the light distribution. Their formalism (the 
so-called IMCAT^) is to date one of the few taking into 
account smearing and anisotropy effects from an atmo- 
spheric PSF. Their formulae were derived for the weak- 
shear limit and assuming that PSFs can be written as 
an isotropic PSF convolved with a compact, anisotropic 
kernel. In this paper we investigate the accuracy and lim- 
itations of this method with simulations. We address the 
following questions: (1) It was shown in earlier publica- 
tions that PSF profiles, for which the KSB breaks down, 
can be constructed very easily (Kuijken 1999). So, is the 
proposed PSF correction procedure valid for PSF profiles 
that we observe in ground-based observations? (2) What 
influence does pixel noise have on the shear estimation? 
(3) The KSB formulae are valid to first order in the shear 
(k <C 1, I7I <C 1). Is this expansion vahd for typical weak 
lensing applications with ground-based observations (from 
I7I = 0.01 for cosmic shear up to I7I = 0.2 for cluster 
mass reconstructions) or should we expand the formal- 
ism to higher order? (4) Can we set up a fully automatic 
procedure, from reduced images to an object catalog, for 
reliable shear measurements? 

We generate a large number of simulated images and 
analyse them exactly in the same way as real data. We 
perform two kinds of simulations: 

— Simulations where all objects involved (galaxies as well 
as stars) have Gaussian profiles. 

— Simulations with the SkyMakcr program (Bertin 2000, 
in preparation) . This program generates galaxies mod- 
elled as exponential disks with a central de Vau- 
couleurs type bulge of varying ratio. The most impor- 
tant feature of SkyMaker is its ability to produce real- 
istic ground-based PSFs and to allow for the inclusion 
of several telescope/detector defects (like astigmatism, 
coma, drift) 

The outline of the paper is as follows: Sect. 2 summarises 
the KSB formalism and how we use it to estimate grav- 
itational shear. Sect. 3 gives a short overview over the 
Stuff program which generates the input galaxy and star 
catalogs for the SkyMaker program. We describe our sim- 
ulations with Gaussian profiles in Sect. 4, our object de- 
tection and selection procedure in Sect. 5. The SkyMaker 
simulations follow in Sect. 6 and we finish with our conclu- 
sions in Sect. 7. The generation of Point Spread Functions 
in the SkyMaker program is described in an Appendix. 

2. Shear estimates 

KSB estimates the shear by considering the first order 
effects of a gravitational shear and an instrumental Point 
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Spread Function (PSF) on the (complex) ellipticity of the 

galaxies X = Xi + 1X2, which is defined as: 



Qu — Q22 



2Q 



Qu + Q22 Qu + Q22 

where the Qij are weighted second brightness moments of 
the light distribution I: 



d^9e,9ji{e)wi\e\), 



(2) 



where the object centre is at the origin of the coordinate 
system, i.e. 



d^eei{0)w{\0\) = 0. 



(3) 



The total response of the galaxy ellipticity x that is the 
intrinsic ellipticity convolved with an isotropic function 
[see Bartelmann & Schneider (2000)] , to a reduced gravi- 
tational shear g = 7/(1 — k) and the PSF is given by: 
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pg _ psh. _ psm ^p*sm'^ — l p*sh 

where x is the observed ellipticity and the tensors P^^ 
and P**™ can directly be calculated from the galaxy's light 
profile and the weight function W. P^^ {Shear polar iz- 
ability) would be the response of the weighted galaxy el- 
lipticity to a gravitational shear in the absence of PSF 
effects. modifies this tensor by a factor including the 
Smear polarizability tensor P*^™ to calibrate the shear es- 
timate for the circular smearing by the PSF. This cali- 
bration also depends on the corresponding tensors P*^^ 
and P***™ from stellar objects containing the information 
of the PSF. The stellar anisotropy kernel q*, needed for 

the correction of the PSF anisotropy, can be estimated by 

A. * 

noting that x^ = 0, g* = for stars so that 



q* _ ^p*sm^-l^* 



(5) 



A complete derivation of these formulae clarifying all the 
assumptions in the formalism is published in Bartelmann 
& Schneider (2000) but see also KSB, LK97 and Hoekstra 
ct al. (1998). To estimate g we first correct objects for the 

PSF anisotropy: 



(6) 



and then consider averages over galaxy images in areas 
where we assume constant reduced shear: 

where we used (x*^) = 0. For simplicity eq. (4) is often 
used in the form: 

g = (P^)-'(X""'^° - XO) ^ (g) = ((P^)-'x^"'^°) (8) 

where we have to assume ((P^)~^X°) = 0- We checked 
that with our procedure described below this assumption 
holds and the estimators from eq. (7) and eq. (8) are in 
very good agreement, although in practice the latter is 
easier to handle. 
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2.1. Practical application of the KSB formulae 

Although the general KSB procedure seems uniquely de- 
termined, there are several possibilities for its practical 
implementation. This has been the subject of discussions 
in recent works (Hoekstra et al. 1998, Kuijken 1999, Kaiser 
1999) on how the KSB formulae should be applied opti- 
mally. One issue is the window function with which 
galaxies are weighted. KSB have chosen a Gaussian with a 
scale tJ proportional to the size of the object under consid- 
eration. Hoekstra et al. (1998) noted first that quantities 
from stellar objects should be calculated with the same 
scale as the objc;ct to be corrected (this is the most natu- 
ral choice although according to the derivation of the KSB 
formulation the result should not depend on it). Another 
difficulty arises as the quantities q* and (p*sm-j-ip*sh 
have to be known at the positions of the galaxies to be 
corrected. Since these quantities are known only at star 
locations, we need to estimate what the PSF would bo on 
the position of each detected galaxy. On top of all, all mea- 
sured quantities are affected by noise in an unknown way. 
Two techniques have been proposed to suppress this noise 
in the shear estimate. First, the noise in the scaling factor 

can be minimised by a fit or by calculating means or 
medians in bins in the parameter space 7rpg = (object size, 
magnitude) which has often been used in the literature. 
Second, we can easily introduce a weighting scheme in the 
estimators (7) and (8). 

Our procedure, which we found to be optimal in terms 
of shear recovery accuracy, can be summarised as fol- 
lows: To apply (4), for every galaxy the quantities q* and 
j^p*sm-j-ip*sh j^g^yQ known everywhere on the image. 

Usually the PSF anisotropy varies smoothly over the field 
so that the q* can be well represented by a low-order poly- 
nomial fit from the light profiles of bright, unsaturated 
foreground stars. A second-order polynomial is enough, 
and higher order polynomials do not improve the PSF 
correction. We select these stars in a rh vs. m diagram 
(see Fig. 1), and then measure q* with a filter of the size 
of the stars, but (p*sm)-ip*sh ^ range of filter scales 
spanning the sizes of our galaxies. We then perform the 
anisotropy correction (6) for each individual galaxy. For 
Ps we finally use the raw, unsmoothed values as our anal- 
ysis shows that fitting this quantity does not improve the 
final results (see Sect. 4). 

Corrected galaxy ellipticities are subject to high noise, 
which can eventually produce unphysical ellipticities much 
larger than unity. It is therefore necessary to weight each 
galaxy according to the accuracy of its ellipticity measure- 
ment / correction. For our weighting scheme we assume the 
ellipticity distribution of the lensed galaxies is to first or- 
der equivalent to the intrinsic ellipticity distribution ev- 
erywhere on our data field. We also assume that the in- 
trinsic ellipticity distribution does not vary over the field. 
So, high ellipticities are assumed to originate from noise 
and these galaxies should be assigned a low weight. These 



assumptions are valid for studies of empty fields but not 

when investigating cluster lenses where high ellipticities 
are lensing features and not necessarily caused by noise. 
In this case a weighting scheme calculating weights out of 
the pixel noise properties may be adopted. See the Ap- 
pendix of Hoekstra et al. (1999) for an example. 

We introduce a simple U = l/cXg weighting, where (Tg 
is the variance of the shear estimators from individual ob- 
jects. These weights are estimated in a parameter space ttu 
supposed to trace the noise properties of objects (like ob- 
ject size and S/N for example). Unfortunately our objects 
are typically very clustered in i^u so that calculating aver- 
ages in cells defined by a regular grid is not optimal and an 
adaptive grid should be more appropriate. Therefore for 
each galaxy we consider its N nearest neighbours in our 
parameter space (typically N ~ 20) and assign the inverse 
of (Jg = X] IsP from these neighbours as weight to that 
galaxy. The distance d from a galaxy k to its neighbour I 
in a parameter space tt consisting of i = 1..M elements is 
hereby defined by 



\ 



M 



(9) 



Obviously, a galaxy does not have an unique set of clos- 
est neighbours in an arbitrary space. For instance, object 
sizes would simply scale by a factor of 2 if data are re- 
binned to half the original resolution while quantities like 
m stay unchanged. Possible consequences of this are not 
investigated here. Hereafter we denote averages (.x) and 
uncertainties of means oi^^) from a weighted quantity x 
by: 



{x) 



•{x} 



1 



(10) 



As noted above, our weighting scheme is well justified only 
if the lensed ellipticity distribution is to first order equiv- 
alent to the unlensed one. We found in our simulations 
that we can also estimate shear that is not small with re- 
spect to the intrinsic ellipticity distribution with the above 
formulae but not the corresponding errors. The reason is 
that only relative weights enter the shear estimate while 
the errors cr^^^ are determined by the absolute values of 
the weights. In the strong shear case, that is. if the in- 
trinsic ellipticity distribution is substantially changed by 
the input shear, the absolute values of the weights directly 
depend on the input shear. Hence we first use eq. (10) to 
estimate the reduced shear according to (8): 



9c. = {{p^Uxfn 



(11) 



and then estimate the errors by bootstrapping. In parallel 
to using the full tensor expression for Ps we also calculate 
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Fig. 1. Plotted is the half light radius vs. the magnitude (in arbitrary units) from all the detections in one of our SkyMakcr 
simulations in Sect. 6. Bright, unsaturated stars are clearly localised in a small region around r?i ~ 2 and in 15 < m < 20. Larger 
objects are considered as galaxies, while the smaller, very faint objects are probably noise detections. For objects around the 
stellar locus with faint magnitudes, a clear classification between stars and galaxies cannot be done and we typically exclude 
this region from our analysis. 



the often used scalar correction, where is estimated by 



Pf = 0.5 tr[PS] ; (g) 



PI 



(12) 



(Hudson ot al. 1998, Hockstra et al. 1998). The weights 
for the tensor and scalar case are generally different. 

3. Catalog generation 

For our simulations we used the Stuff program^ to gener- 
ate our initial object catalogs. Details of this program and 
the SkyMaker tool producing images out of these catalogs 
will be published elsewhere (Bertin 2000, in preparation); 
here, only a very short summary of those parts relevant 
for this work are given; 

Given the sky dimensions of the intended simulations, 
the program distributes galaxies in redshift space that is 
subdivided into bins. For every bin, representing a vol- 
ume element according to the specified cosmology (we 
used Ho = 65km/(s Mpc), fim = 0.3 and fix = 0.7 in 
all our simulations), the number of galaxies for the Hubble 
types E, SO, Sab, Sbc, Scd and Sdm/Irr is determined from 



Freely available 
ftp : //ftp . iap . fr/pub/f rom_users/bert in/stuff / 



at: 



a Poisson distribution assuming a non-evolving Schechter 
luminosity function (Schechter 1976). The different galaxy 
types are simulated by linearly adding exponential iJ,d{r) 
(disk component) and de Vaucouleur profiles Hh{r) (bulge 
components) in different ratios. 



1/4 

^b(r) = Mb + 8.3268 ( — ) + 51ogrb + 16.6337 (13) 
/Zd(r) = Md + 1.8222 ^-^ ) + 5 log ra + 0.8710, (14) 



where fXh, Hd are the surface brightnesses in mag/pc^, and 
Afb, Md are the absolute magnitudes of the bulge and the 
disk components, respectively. The distributions of scale 
radii rb and ra are fixed by an empirical relation (Binggeli 
et al. 1994) and a semi-analytical model (de Jong & Lacey 
1999) relating these quantities to the absolute magnitudes. 
The galaxies are assigned a random disk inclination an- 
gle and a position angle which define the intrinsic ellip- 
ticities of our objects. The output of the program is a 
catalog of galaxy positions, apparent magnitudes, semi- 
minor and major axes, and position angles for disks and 
bulges. These catalogs are then sheared (disks and bulges 
are sheared separately!), processed with SkyMaker and 
run through our KSB procedure. Before describing these 



How accurately can we measure weak gravitational shear? 



5 



SkyMaker simulations we present results with Gaussian 
object profiles in the next section. 

4. Semi-analytical calculations with Gaussian pro- 
files 

Although Gaussian profiles do not fit real galaxies and 
stars very well, they allow a quick investigation of the 
possible biases connected with the KSB procedure that 
arise from constraints in real data. Moreover most of the 
quantities defined in KSB can be analytically calculated 
with Gaussian profiles, which allows one to check that the 
numerical simulations are done properly and to explore 
the following effects: 

— CCD images are pixelised and all integrals involved in 
the KSB procedure have to be estimated by discrete 
sums. 

— Our data are affected by sky noise. The expressions 
from Sect. 2 are non-linear in the object profile 7, so 
noise in this quantity may introduce systematic biases. 
Sky noise also affects the calculation of the objects 
centroid position from which all shear quantities are 
calculated. So wc have to investigate how the noise in 
object positions influences the final shear estimates. 

— In practice we only have the observed light distribu- 
tion for our analysis. When applying the KSB formal- 
ism we make the assumption that the influence of the 
anisotropic component of the PSF and shearing on the 
profiles arc "small" , although in principle intrinsic pro- 
files should be used. 

Furthermore we investigate different possible parameter 
spaces TTpg and wu in view of a possible noise reduction 

in Ps and our weighting scheme. Our Gaussian galaxy 
profiles were convolved with Gaussian PSFs given by: 



/ = Aexp(-- 



62 



(15) 



where x' = xcos{9) — ysm{9), y' = xsm{6) + ycos{9) for 
an object with semi-major (semi-minor) axes a and b and 
a position angle 9 with respect to the abscissa. The dimen- 
sions of the pixels and the PSF profiles were adapted to 
a CCD with 0'.'2 resolution and a seeing (FWHM) of Cf.'l 
(This means a* = 5* = 2.012 pixel imits for an isotropic 
PSF). We chose A* = 10000 to get high S/N measure- 
ments for all stellar quantities. The intrinsic scales for the 
galaxies a and b were taken from the output of the Stuff 
program. The distributions of these scale radii are very 
peaked at around 0.2" and have a long tail towards larger 
radii. The formal mean intrinsic size (defined as \^ab) of 
the galaxies is 0.22" ± 0.13". The apparent magnitudes 
from Stuff were used to fix the amplitude A. The centres 
of the objects were chosen randomly within a pixel. 

We then performed the KSB measurements on these 
profiles for a noise-free, and for a Gaussian noise model 
with (Tsky = 13.0 (this adapts these calculations to the 



characteristics of the simulations in Sect. 6.1). Hereby 

we added the noise only to the galaxies but not to the 
stars. So we can investigate the best possible noise-free 
results only having the pixelisation as a source of error, 
and the consequences of sky noise on them. In the follow- 
ing analysis we used all the galaxies in the initial Stuff 
catalog regardless of whether the objects would be found 
by some detection algorithm in a real image or not. For 
the weight function Vl^(|6'|) we chose here a Gaussian with 
scale Tg = 1.107 , where as and 6s are the scale 

radii of our objects after smoothing with the PSF. This 

corresponds to the rg radius from KSB (see Sect. 5) for 
a circular Gaussian as = 6s. We have checked that the 
results that we present in this section do not depend sig- 
nificantly on this choice. Before calculating any lensing 
quantities we estimated the object centre by iteratively 
solving the equation: 



j9I{9)W{\0-9o\)d'^9 

J i{e)w{\e - eo\)d^e ' 



(16) 



where we used the pixel centre of the true position as start- 
ing point. In addition, for every object we have defined a 
signal-to-noise ratio by: 



_S _ J I(9)\V{\0\)d'^9 

N <7sky^/W^2(|^|)d2^' 



(17) 



as we wanted a S/N ratio that is based on the same filter 
function as the one that the measurements are done with 
(in the end it turned out that there is a linear relation 
between our S/N estimator and the v estimator from the 
KSB hierarchical peak finder; see Sect. 5). For the sample 
with S/N> 2 the rms of Ax and Ay in the determination 
of the object centre is better than 1/20 pixel. For the low 
S/N< 2 objects the rms is about 1 pixel in both com- 
ponents. All other quantities were then calculated with 
respect to the estimated centre, where the integration was 
done up to |Aa;| < 3rg; \Ay\ < 3rg. We generated cat- 
alogs for the seven shear combinations listed in Table 1 
that cover the range from |g| = 0.012 to |g| = 0.32. For 
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0.01 


0.03 


0.15 


0.25 


0.08 


0.2 


0.14 


92 


0.007 


0.05 


0.1 


0.2 


0.27 


0.14 


0.2 



Table 1. The shear combinations we investigated in this sec- 
tion with Gaussian profiles for the galaxies. 



each combination wc generated a final catalog with about 
75000 objects having 15.0 < m < 28.0. If not stated oth- 
erwise our galaxies have the intrinsic ellipticity generated 
by the Stuff program. 
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4.1. Estimation of 

As we have calculated all quantities with and without sky 
noise we can investigate whether the smoothed in some 
parameter space ttps can lower the noise in this quantity, 
and whether this improves the final shear estimation. For 
this we took the gi = 0.25; , 92 = 0.2 realisation and the 
scalar PS estimator from eq. (12). We smoothed P^ with 
the next neighbour approach described in Sect. 2.1 in sev- 
eral parameter spaces to obtain -P^„ooih and compared 
these values with the noise free P^^^^. Here we finally used 
the median to estimate -Pgniooth ^^^^ turned out to give 
better results than the mean. The result is shown in Fig. 
2. We see that P^ is mainly a function of the object size fg 
and the modulus of the raw, noise-free but unobservable 
ellipticity Xtruc- In practice, we only observe noisy values 
for X and a smoothing as a function of rg alone seems as 
good as smoothing in rg and |x|- Including quantities like 
m or S/N is clearly worse. 

In order to see whether these smoothings improve the 
final shear estimates, we first calculated the estimators 
g = -p^ and g^mooth _ — ^^j. ^j^g various smoothings. 

smooth 

Fig. 3 shows that even for the best rg — |xtruc| wc find 
Ag « ^(^smooth^ This shows that smoothing docs not im- 
prove the shear estimates on the whole but we still have 
to investigate whether it helps for objects having an in- 
trinsically unphysical negative P^. For this we compared 
{Ag) and ^A(/™°°*'') together with the corresponding un- 
certainties for the subsamples with P^ > and P^ < 0. 
To exclude extreme outliers from this calculation we in- 
clude only objects with S/N> 2.0 and a modulus smaller 
than unity for the corresponding shear estimate. We ob- 
tain (A(/'"°°^'^) = 0.028 ± 0.19, (Ag) = 0.028 ± 0.18 
for the PS > sample, and (Ag™°oth) ^ 0.16 ± 0.38, 
(Ag) = 0.26 ± 0.42 for Ps < objects. This confirms that 
smoothing does not at all change the estimate for P^ > 0. 
It also shows that we cannot recover the shear signal with 
ps < objects, whether we use a smoothed value for this 
quantity or not. 

4-2. Weighting scheme 

Similarly, we now search for the best possible parame- 
ter space to calculate cTg, and thus the weighting of the 
galaxy ellipticities. For this investigation we used a catalog 
with intrinsically round objects to isolate pixel noise er- 
rors from those of the intrinsic ellipticity distribution. The 
input shear is gi = 0.03, 52 = 0.0. Fig. 4 shows that us- 
ing TTu = (rg, m) and ttu = (vg, S/N) as parameter spaces 
for the (j^ calculation give the best and very comparable 
results in the final shear estimation. The errorbars in this 
and all following calculations are estimated by creating 
200 bootstrap realisations of our catalog in the follow- 
ing way: We randomly take N objects out of our catalog 
containing N objects allowing for multiple entries. Prom 
this sample we recalculate a shear estimate and repeat the 



exercise A/boot = 200 times. From these shear estimates 

^boot.^ G [l-.Aboot] the errorbar for the original measure- 
ment' is estimated by a^^^ = ^^.(gboo* - {3"°°°'))^ 

(we checked that {gf°°^) is in excellent agreement with 
our original shear measurcmcnit). We note from Fig. 4 that 
introducing our weighting scheme improves the shear esti- 
mation significantly in comparison to no weighting. There- 
fore, weights from the ttu = (rg,S/N) smoothing are now 
used throughout. 

4.3. Accuracy of shear estimates 

To evaluate how accurately we can measure gravitational 
shear we have investigated the gi, 92 combinations listed 
in Table 1. We did our measurements first without PSF 
anisotropy and then with a stellar axis ratio 6* /a* = 0.9 
with a* oriented such that (73 = 0. To obtain the results 
quoted below, we have calculated averages (g) and er- 
rors cr(g) over all galaxies in our catalogs. The results for 
the noise- free calculations are shown in Fig. 5. Note here 
that we have also used the same weighting scheme for the 
noise-free calculations, although this would not be neces- 
sary. From the plot we see the following trends for the 
calculations without PSF anisotropy: 

— For small input shear values |g| < 0.1 we can clearly 
recover the input shear within Ag « 0.01 with the 
tensor as well as with the scalar shear estimator. For 
larger shear values we systematically underestimate 
the shear, with a maximal underestimation of 0.02 to 
0.03 for the highest input shear of 0.2 — 0.25. We note 
that for both shear estimators the underestimation as 
a function of input shear can be very well represented 
by a straight line indicating that the underestimation 
is a constant fraction of the input shear over the whole 
range of input shears considered here. This fraction is 
about 10% - 15% for the scalar and 5% - 10% for the 
tensor shear estimate. 

— We also note that the results are completely equivalent 
in both shear components. 

For the simulations with PSF anisotropy in the gi direction 
we conclude the following: 

— Considering the g2 component the results are very sim- 
ilar to the isotropic PSF case. 

— For the gi component the scalar correction underes- 
timates the shear by an amount comparable to the 
isotropic PSF case. In contrast, the tensor correction 
overestimates the input shear in the direction of the 
PSF anisotropy. This relative overestimate is about 
20% for the high |g| and can reach up to 100% for the 
lowest input shears. In Fig. 8 we show analytic results 
with intrinsically round objects and an input shear of 
gi = 0.25. These calculations confirm that for Gaus- 
sian profiles the tensor correction overestimates the fi- 
nal shear in the direction of a PSF anisotropy. From 
our analysis of the SkyMaker simulations in Sect. 6, 
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where we have encountered similar problems with the 
tensor correction, we conclude that the problems are 
not connected with the anisotropy correction of KSB 

but with the smearing correction. 

For the calculations including sky noise we have investi- 
gated the whole sample, and two sub-samples with S/N< 
2 and S/N> 2. Fig. 6 shows that the S/N> 2 sample 
results are similar to the noise-free case, i.e. in the ab- 
sence of PSF anisotropy we have very good results with 
the tensor estimate (however, there is no underestimation 
in the tensor case at all in contrast to the noise-free case; 
in fact, the underestimation is lowered to 4% in the noise- 
free case for the same subsample of galaxies), and about 
a 10% — 15% underestimate for the scalar calculations. 
For the low signal-to- noise objects (S/N< 2) the shear 
is underestimated by about 30% in the scalar and up to 
about 50% in the tensor case, where the underestimate in 
the tensor case shows an increasing trend for increasing 
shear. The error in the object centre determination is not 
the main reason for the underestimate: for S/N< 2 we re- 
peat the same calculations using the true object centres. 
The results from this test arc very comparable to those 
shown here. For the whole sample, the scalar correction is 
better with a relative underestimate of about 10% — 15% 
over the tensor one which has an underestimate of about 



20% — 30%. If PSF anisotropy is present, noise-free and 
noisy cases give the same trends. Fig. 7 shows the noise 
amplitude of the estimators. We considered the ratio of 
the shear uncertainty from noisy objects a and the noise 
free calculations a^^^^. For S/N> 2 pixel noise increases 
the errors in the shear estimates by about a factor of 2, 
for S/N< 2, by a factor 5-10, and by a factor of 5 for the 
whole sample. 

4-4 ■ Summary 

With our proposed way to apply the KSB technique to 
measure gravitational shear we come to the following con- 
clusions when dealing with Gaussian profiles: 

— Our analysis shows that the scalar estimator from 
eq. (8) systematically underestimates shear by about 
10% — 15%. For low-S/N objects, the underestimate is 
about 30%. The tensor estimator from eq. (7) is better 
for high-S/N objects showing no systematic over- or 
underestimation in the final result. This estimator is 
less stable for low-S/N objects where there is a system- 
atic imderestimation that can reach more than 50%. 
When including PSF anisotropy, the scalar estimator 
still gives very stable and comparable results both in 
the direction of the PSF anisotropy and perpendic- 
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smoothing of this quantity with rg and the true raw cUipticity Xtruc we have Ag ~ ^^smooth jjjfjjgating that smoothing 
does not improve the shear estimation. The parameter spaces rg and rg — \x\ give similar results while the rg — S/N smoothing 
especially shows a spread of /\gf™°°*'' where Agi = 0. The result implies that smoothing does not improve the final shear 
estimate. We checked that this general result also applies for objects with an intrinsic, unphysical negative value for P^ (see 
text). 



ular to it. In contrast, we overestimate the shear in 
the direction of PSF anisotropy in the tensor case. We 
conclude that the scalar estimator is more stable and 
conservative. Both estimators show very similar noise 
properties. 

— We have shown that smoothing does not improve 
the final shear estimates over taking raw, noisy values. 

— We can give an objective parameter cut S/N> 2 for 
which we can measure shear with about the same ac- 
curacy as with no sky noise. Our S/N parameter [see 
eq.(17)] is very convenient as we do not have to find the 
threshold for every observation individually (as would 
be the case for quantities like m). 

As a final step we have repeated the analysis done here 
putting our Gaussian galaxies and stars on FITS images 
and analysing them in the same way as the SkyMaker 
simulations described in the next section. In addition to 
the steps presented so far, we first have to perform object 
detection and selection for the shear determination. The 
steps for this are summarised below in Sect. 5. An im- 
portant difference with the previous semi-analytical cal- 
culations is the fraction of S/N< 2 objects. Previously 
they constituted more than 60% of all the objects, but for 
the simulations where objects have to be detected first, 
this fraction is very small, about 5% — 10% in the final 



catalogs. It turned out that shear estimates with the sub- 
sample S/N> 2 and the complete sample are now similar, 
so that the small subsample of low-S/N objects becomes 
unimportant. Therefore, only the results with all objects 
without a cut in S/N will be shown from now on. Fig. 9 
displays the results from the simulations. They are com- 
parable to those obtained from the S/N> 2 subsamples 
in Fig. 6. We also conclude that the stellar objects de- 
tected in the FITS images have sufficicint signal to noise to 
give reliable estimates for P^. As described above we only 
used very high S/N stellar objects in the analytic calcula- 
tions. This last simulation ensures that our shear analysis 
pipeline gives results consistent with analytical prediction, 
and that it is ready to be used for the final analysis in- 
volving the realistic simulations done with SkyMaker. In 
the next section, we first describe our procedure leading 
from the image frames to a galaxy catalog for shear mea- 
surements. Afterwards we will see whether our results are 
still valid when considering different galaxy profiles and 
especially more realistic PSFs. 

5. Object detection and selection 

In this section we describe a fully automatic procedure to 
go from image frames to a final object catalog for shear 
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Fig. 4. The measured reduced shear gi as a function of the object size for various parameter spaces used in our weighting 
scheme. Panel (a) shows the result for noise free objects. Panels (b) and (c) show results if no weighting at all is applied, (d) 
and (e) for rg and rg— m smoothings and finally (f) and (g) for the parameters rgt;— S/N. The error bars are from 200 bootstrap 
realisations of the galaxy samples (see text). We note that the weighting scheme improves the measurement significantly to 
the cases where no weighting scheme is applied. We also conclude that there is no need to exclude low S/N objects from the 
measurement if weighting is applied. 



analysis. We note that our main intention was to obtain a 
catalog ensuring reliable measurements. We demonstrated 
at the end of the last section that our analysis of isolated 
objects in the semi-analytical treatment and with those 
detected in FITS images gives very comparable results. 
Our procedure contains several conservative rejection cri- 
teria for objects and it is not optimised to make the max- 
imal use of data in terms of number density of objects. 
The procedure consists of the following steps: 

1. Objects were detected with the 'hfindpeaks' algorithm 
from Nick Kaiser: The image is smoothed with Mexi- 
can top hat filters Wd(|^|) 



1 



2rh 



exp 



2rl 



1 

-o exp 



2r1rf 



(18) 



with increasing filter radii rg. In every smoothed im- 
age, the peaks are detected and linked with the peaks 
found in previous smoothings. Hereby peaks are as- 
sumed to belong to the same object if their positions 
coincide within Tg. In this way a peak trajectory is 
built up for every potential object. For every peak, a 
signal-to-noise ratio v is calculated where the signal is 
the peak value in the smoothed image and the noise 



CTd is given by 



Csky 



= CTsky 




(19) 



with rf = 2. Fig. 11 shows a comparison between this 
signal-to-noise estimate and ours from eq. (17). Af- 
terwards, peaks with the highest i/max value for each 
trajectory identify objects (peaks with fmax < 4 are 
immediately rejected). This size rg is used in all the 
following analyses and the pixel centre of the peak po- 
sition is taken as a starting point for the object centre 
determination. To get a more accurate object centre, 
a Newton-Raphson step is performed and the x and y 
pixel-positions are corrected from the pixel centre by 
5x and 6y 



Sx = - - 



5y = 



(20) 



where J|, and I^^, lyy are the first and second 
derivatives of the smoothed light profiles along x and 

y- 
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direction. No object selection was done and our weighting scheme was used. The points were obtained by averaging over all 
w 75000 galaxies that we have for every shear realisation in our catalog. The error bars are less than 0.1 percent and so smaller 
than the symbols. For clarity of the plot the errorbars have been omitted and the symbols for the scalar correction have been 
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2. From the peaks (objects) found in the previous step, 
all relevant quantities (like P^^, P'^™) are calculated 
separately with the original, unsmoothed image. 

3. From the catalog generated in the previous step, we 
first removed objects with obvious problems during 
the analysis process: (1) Objects closer than 3rg to 
the border of the image; (2) objects where one of the 
eigenvalues A^, B'^ from the Qij tensor was negative. 



P2 




+ 4Q?2;(2i) 



(3) the object has a total negative flux. (4) Fig. 10 
shows the distribution of the difference between the 
centroid estimates of 'hfindpeaks' and ours from eq. 
(16). All objects with a difference larger than 0.1 pixel 
were excluded. We found that this is an efficient way 
to reject blended objects, since one centroid estima- 
tor is based on the smoothed and the other on the 
unsmoothed image. 

4. Objects with a neighbour within 3rg are rejected. 

5. Stars are preselected using the star branch of the rh— m 
diagram and polynomials for the two components of q* 



are calculated in the following way: a preliminary fit is 
done for ql and using a minimisation 



(22) 



AT* 

Xj = ^(.qj{xi,yi)-Pj{ak,Xi,yi)f, 

i=l 

where N* is the number of preselected stars, j = 1,2, 
Xi,yi are the positions of the stars and pj are two- 
dimensional, second order polynomials with six un- 
known parameters afe. After determination of ak we 
calculated the expected error ag^ for every g| by 



CTa* = 



N* 



(23) 



where . is the minimum of at the fitted param- 

-^^J min ^ ^ ^^J 

eters a^. Stars which deviate at more than Icr^^ in any 
of the two components qi or q2 are rejected, and the 
fit is repeated for the final polynomials. 
6. The final sample of stars is reprocessed nine times 
with filter scales = 1..9, in order to match all the 
possible galaxy sizes. After each processing, objects 
with problems according to step (3) are still rejected. 
From the remaining stars we calculated the mean of 
(tr[P*«'^]/tr[P*^"]) and we used this value for the PSF 
correction. 
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Fig. 6. The same as Fig. 5 for noisy objects. The two upper most panels show results for S/N> 2.0, the middle panels for 
S/N< 2.0 and the lower panels for the whole sample. For high S/N objects, the tensor estimation is more accurate but in general 
the scalar calculation is much more stable towards low signal objects. Symbols for the scalar correction are slightly shifted to 
the right for claxity of the plot. For more details see the text. 



7. We now have all the quantities in hand to calculate 
scalar and tensor shear estimates for every objects as 
described in Sect. 2.1. Hereby we considered only ob- 
jects as galaxies whose half light radius was larger than 
the stellar locus (see Fig. 1). 

8. Weights for the tensor and scalar g are calculated, 
which is an important ingredient for the two shear 
estimators as described in Sect. 2.1. With the whole 
procedure we end up with a number density of about 
30 galaxies per sq. arcmin. 



seeing of (/.'? is shown in Fig. 12. A description how the 
PSF is constructed in SkyMaker is given in the Appendix. 

6.2. Results of the SkyMaker simulations 

Fig. 13 shows contours of the outer and core parts of the 
PSFs we have used in our simulations. The outer parts 
of the profiles all look the same (except for the quadratic 
PSF) and they differ mostly in the cores. The anisotropies 
caused by these PSFs are given in Table 2. For every PSF 



6. The SkyMeiker simulations 

6.1. The image characteristics of our SkyMaker simula- 
tions 

For every combination of the shear and PSF, ten Sky- 
Maker images were produced and analysed. Each image 
has the characteristics of a 12000 sec. / band exposure 
taken at a telescope with a 3.5m primary and a Im sec- 
ondary mirror, and has a dimension of 2048 x 2048 pixels 
with a scale of 0'.'206 per pixel. The atmospheric seeing 
in every image was 0'.'7. For every image we detected our 
objects and analysed them as described above. A compar- 
ison of magnitude distribution, object size and ellipticity 
distribution between detected objects in one of the images 
and a 12000 sec. / band CFHT image having a measured 
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PSF_3 


PSF_4 
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0.00 


0.00 


0.00 


X2 
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0.01 
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PSF_7 


PSF_8 


Xi 


0.1 


0.06 


0.04 


0.00 


X2 


0.07 


0.03 


-0.03 


0.00 



Table 2. The anisotropies of the PSFs described in Fig. 13. 

The values quoted are the raw ellipticity measurements from 
KSB where the from the stars was used as filter scale. 



we have analysed seven sets of images with the shear com- 
binations given in Table 3 where a few sets are slightly dif- 
ferent from those in Table 1 . The results for our SkyMaker 
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Fig. 7. The noise properties of the shear estimations from Fig. 6. For high S/N objects, sky noise only has a very weak 
contribution to the shear estimates. It increases the uncertainty by a factor 3-5 for low S/N objects and by a factor of about 
2-3 for the whole sample. Symbols for the scalar correction are slightly shifted to the right for clarity of the plot. The panels 
are for the same object samples as in Fig. 6. 



91 


0.01 


0.03 


0.25 


0.08 


0.2 


0.2 


0.14 


92 


0.007 


0.05 


0.2 
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Table 3. The shear combinations we investigated with our 
SkyMaker simulations. 



simulations arc shown in Figs. 14, 15, 16 and Table 4. As 
for the semi-analytical and simulated Gaussian profiles we 
now discuss the accuracy of the shear estimators for the 
different shear amplitudes and PSF anisotropies: 

— In the case of no PSF anisotropy at all (PSF_1), the 

tensor correction gives the correct shear on the whole 
range of |g| = 0.012 up to |g| = 0.32. The scalar cor- 
rection underestimates the shear relatively by about 
10% - 15%. The la error bars are about at ±0.005. 

— In the presence of PSF anisotropy, the tensor correc- 
tion can over- or underestimate the true shear while 
it is always underestimated in the scalar case. In the 
sense of deviation from the true value, the tensor case 
always gives the better result. The relative underesti- 
mates with the scalar correction can reach up to 30% in 
the worst case. As the anisotropy of stars is perfectly 
corrected by our polynomials for p, we investigated 
whether there is nevertheless an anisotropy residual in 



the galaxies. For this we compared the position angle 
from the input shear with the recovered one in Fig. 
16. There we see that in the scalar case the position 
angle is nearly perfectly recovered. This means that 
the PSF anisotropy is corrected very well but the cal- 
ibration factor Ps is too low. In the tensor case the 
recovery of the position angle is not as good as in the 
scalar case but also very acceptable. The reason for 
the over- or underestimates of the magnitude of the 
shear is probably that the elements of that are used 
for the isotropic correction should be calculated from 
profiles that are already corrected for PSF anisotropy 
defects. 

— From Table 4 we conclude that also for all the Sky- 
Maker simulations, the over- or underestimates are 
nearly a constant fraction of the input shear. For the 
tensor estimator we can typically recover the input 
shear to an accuracy between 10% — 15%:. This means 
that we can measure a shear below 0.1 nearly within 
0.01, the accuracy required for accurate measurements 
of the cosmic shear on scales < lO'.O. 

— For all the SkyMaker simulations the recovery of the 
shear is better than in the anisotropic case with Gaus- 
sian profiles (see Fig. 9). It seems that the more real- 
istic SkyMaker PSF profiles better reflect the assump- 
tions in the KSB algorithm than Gaussian profiles, 
namely that PSF anisotropy mainly comes from the 
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Fig. 8. Analytical results for gj = 0.25, round intrinsic objects, no PSF anisotropy and a 10% PSF anisotropy in the gi direction 
are shown. The solid and short dashed lines are scalar and tensor shear estimate for the PSF anisotropy free case, the long 
dashed and dashed dotted line represent scalar and tensor correction when PSF anisotropy is present. The figure confirms that 
the resulting shear is overestimated significantly in the tensor case with PSF anisotropy. 



central core. While a Gaussian has its anisotropy on 
all scales, Fig. 13 shows that the anisotropy details of 
the SkyMaker PSFs lie in the core indeed. 
— Regarding the issue of scalar vs. tensor correction, we 
can conclude that, although the scalar case always un- 
derestimates the true shear, it provides the more con- 
servative answer and it should be used when one is in- 
terested mainly in the position angle of objects while 
the tensor should given preference when the amplitude 
of the shear is important. 

7. Conclusions and Outlook 

In this paper we have investigated how well we can re- 
cover weak gravitational shear with realistic simulations 
for ground-based observations. With respect to the main 
motivation for this work, to test the reliability of our re- 
cent detection of cosmic shear, the results are very en- 
couraging. In all our simulations we could recover weak 
shear up to |g| = 0.1 with an accuracy of 0.01 or bet- 
ter and we can significantly exclude the detection of a 
I7I = 0.04 — 0.05 false signal based purely on uncorrected 
PSF effects. Although the PSFs we have investigated do 
certainly not cover all possible defects mimicking a lens 
signal, also other independent studies like Hoekstra et 
al. (1998), who looked at diffraction- limited model PSFs 
for the HST or Bacon et al. (2000b), indicate that the 



KSB algorithm works better than one could expect look- 
ing at its assumptions. We showed that for the PSFs in- 
vestigated here, the anisotropy correction proposed in the 
original KSB work is working very well, but that there 
are problems with the isotropy correction if an anisotropy 
is present. The reason for this is probably that in such a 
case, the boost factor is calculated with the wrong, not 
anisotropy-corrected, galaxy profile. Moreover, we have 
presented a fully automatic procedure leading from the 
image frames to an object catalog for reliable shear mea- 
surements. Despite these encouraging results, the advent 
of new wide-field imaging mosaic cameras in particular, 
and the associated data flow brings additional technical 
problems if we want to build up a fully automatic data 
processing pipeline. Some of these problems that we could 
not address in the scope of this work, but leave to a future 
publication are: 

1. In wide- field imaging especially, many images suffer 
from very bright stars causing blooming effects, strong 
stray light and reflections. So far, we usually have 
marked out affected image regions by hand, but it 
should be investigated whether there are automatic 
ways of dealing with it. 

2. We usually coadd several single exposures to get a fi- 
nal image. These frames typically have slightly differ- 
ent seeing disks. In the past we often simply did not 
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Fig. 9. Results for Gaussian profiles put on FITS images and analysed with the procedure outlined in Sect. 5 are shown. Panels 
(a) and (b) are for the PSF anisotropy free case and the panels (c) and (d) for a PSF anisotropy of 10% in the gi direction. 
All galaxies were used and no selection according to S/N was done. Dots with errorbars represent measurements with tensor 
estimation, lines with errorbars those with scalar estimation. The final estimates were done with about 13500 galaxies for every 
realisation and error bars have a typical size of 1%. The results are completely comparable to the S/N> 2 case of Fig. 6 



use images that had the worst seeing. We have to in- 
vestigate what effects the coaddition of images with 
different PSF properties has on the final shear mea- 
surement result. 
3. As a consequence of the last point, weak-lensing ob- 
servations have mostly been done with a very compact 
dither pattern so far. This has the advantage of al- 
lowing a simple coaddition of the images with integer 
pixel shifts, and the effects of optical image distortions 
do not need to be taken into account during the image 
prereduction. This first preserves uncorrelated noise in 
the image pixels, second it preserves the smoothness of 
PSF anisotropy on the scale of single chips, and third 
it brings the practical advantage that we can still deal 
with single chips instead of very large images. Never- 
theless, we think that this approach is not ideal. On 
the one hand accurate astrometry and photometry is 
done most easily if information from objects in the 
overlap between different chips is at hand. The afore 
mentioned advantages may turn into problems in later 
analysis when the gaps between single chips lead to se- 
vere border effects, e.g. when we try to search for fila- 
ments between galaxy clusters (see Kaiser et al. 1999). 
So we will investigate whether we can overcome the 
effects of correlated pixel noise caused by remapping 



for optical distortions and discontinuities of the PSF 
anisotropy in the final coadded images. 
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Fig. 10. The distribution of the positional diflferences of the 'hfindpeaks' position estimator and ours from eq. (16). We see a 
clear peak below 0.1 pixel and so we conservatively excluded all objects with a higher difference. See the text for more details. 
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line and our criterion for a good object (S/N> 2) corresponds to f w 7. 
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Fig. 12. A comparison of galaxy properties in one of our SkyMaker images with real data in one of our I band CFHT fields. 
Both images had a seeing of 0'.'7 and both corresponded to 12()()() see. exposures. Objects have been detected and the quantities 
have been calculated with SExtractor. The two upper panels show a comparison of the magnitude distribution, the two middle 
panels of the object size defined as \/ AB, where A and B are defined as in eq. (21) and the lower panels of the ellipticity 
parameter e = (1 — B/A)/{1 + B/A). The only noticeable strong difference is that the ellipticity distribution in the real data is 
slightly broader. The shift of the peak in this distribution is caused by a strong PSF anisotropy in the CFHT data. 
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Fig. 13. The figure shows the outer part (left) and the core (right) of the PSFs used in the SkyMaker simulations. All the 
PSFs contain four spider arms produced by the support of the secondary mirror. Prom left to right and from bottom to top the 
PSFs contain the defects (For the tracking errors, the first value in the bracket gives the rms drift in the direction of the second 
value. The third value is the rms of a Jitter in the orthogonal direction. For the aberrations the "d80" diameters introduced are 
quoted (see the appendix for a detailed description of tracking errors and aberrations). We note that the "d80" diameters of the 
individual aberrations do not add quadratically as the aberration polynomials are not orthogonal.): (1) None; (2) tracking error 
(0.0", 97.0°, 0.05"), defocusing (0.05"), coma (0.26"); (3) tracking error (0.05", 0.0°, 0.01"), defocusing (0.07") and triangular 
aberration (0.015"); (4) tracking error (0.02", 34.0°, 0.02"), defocusing (0.1") and spherical aberration (0.3"); (5) tracking error 
(0.3", 31.0°, 0.2"), defocusing (0.3") and astigmatism (0.72"); (6) tracking error (0.3", 62.0°, 0.3"), defocusing (0.3"), coma 
(0.31") and triangular aberration (0.72"); (7) tracking error (0.5", -20.0°, 0.11"), defocusing (0.07") and triangular aberration 
(0.122") (8) defocusing (0.05") and quadratic aberration (0.61") 
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Fig. 14. From top to bottom the figure shows SkyMaker results for the first four PSFs of Fig. 13. The objects from the ten 
images for every shear/PSF combination were pasted together and the means and errorbars shown here were obtained from this 
pasted catalog. 
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Fig. 15. The same as Fig. 14 for the second four PSFs of Fig. 13. 
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Fig. 16. The figure shows the difference — 0'j of the recovered position angles vs. the position angles of the input shear. 

Crosses correspond to the results from Fig. 14 and diamonds to those from Fig. 15. Triangles are for shear values where no PSF 
anisotropy correction has been done (PSF_4). In this plot we omitted the very weak shear case (pi = 0.01, g2 = 0.007) as the 
position angle is badly determined here. 
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Table 4. The table gives a summary of the results from the SkyMaker simulations performed. Given is the relative error (51,2 — S'i^2°)/5i^2° in percent. We note that 
this error is always a fairly constant fraction over the whole range of input shears (excluding the very weak shear case). The mean number of used objects per sq. 
arcmin and the measured stellar rh are from the catalogs of the gi = 0.01, g2 = 0.007 input shear. 
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Appendix: The SkyMaker program 

SkyMaker is an image simulation program, a kind of "vir- 
tual telescope" originally designed to assess SExtractor de- 
tection and measurement performances (Bertin & Arnouts 
1996). The code (Version 2.3.3^) has been much improved 
since, and is currently capable of simulating star and 
galaxy images with higher accuracy. 

The optical Point Spread Function (PSF) 

Simulated images feature a PSF which is assumed to be 

the convolution of the telescopic (instrumental) point- 
spread function, tracking errors, and atmospheric "see- 
ing" (assuming a long exposure time). We shall however 
more conveniently deal with the Optical Transfer Function 
(OTF), which is just the Fourier transform of the PSF. 

Telescope 

Although the diffraction spot of a large telescope is small 

compared to the atmospheric seeing disk at optical wave- 
lengths, features such as diffraction spikes (created by the 
spider-arms of the secondary mirror /primary focus) are 
obvious on real images and must be taken into account. 

The principle of the PSF generation is similar to that 
of other telescope image simulators like TIM (Hasan & 
Burrows 1995). The system is assumed to be illuminated 
with incoherent, quasi-monochromatic planar waves. In 
those conditions, the instrumental OTF is proportional to 
the autocorrelation of the input pupil function of the sys- 
tem (e.g. Born & Wolf 1999). The on-axis pupil function 
P{p, 0) used in this work is shown in Fig. 17. It is typical 

^ Freely available at 

ftp : / / ftp . iap . f r/ pub/ f rom_users/bertin/skymaker/ 



of telescopes like the 3.6m CFHT, harbouring wide-field 

instruments used for weak- lensing experiments. 

Instead of using ray-tracing techniques to simulate off- 
axis aberrations, we preferred to approximate the latter by 
modulating the phase (f) of the complex on-axis input pupil 
(as shown in Fig. 17). Aberrations are limited to the sum 
of low-order Seidel (1856) terms: defocus 0d cx p'^; astig- 
matism (j)a oc p'^cos^{6 — da); coma (pc p"^ cos{9 — 9c); 
spherical (ps fx p*; plus the triangular and quadratic aber- 
rations. Each term is normalised in "d80" units, that is, 
the diameter of the disk within which 80% of the light from 
a point source image is enclosed, omitting the contribu- 
tion from diffraction and other aberration terms (Fouque 
and Moliton 1996, private communication). 

Tracking errors 

A tracking drift of angle d in the 6t direction can be simu- 
lated by multiplying the OTF by sine {df).d{0 — 9T), where 
/ is the angular frequency. "Jittering" effects are easily 
added by multiplying the OTF by a centered, 2D Gaus- 
sian window. 

Atmospheric turbulence 

Close to zenith, atmospheric turbulence blurs long- 
exposure images in an isotropic way. Experimentally, the 
"atmospheric" OTF is found to quite accurately follow 

the theoretical prediction cx exp —3.442 (j^^ for a 
turbulent atmosphere that follows Kolmogorov statistics 
(see for instance Roddier 1981 and references therein, or 
more recently Racine 1996). Fried's (1966) ro parame- 
ter can be more conveniently expressed as a function of 
the "atmospheric" Full- Width at Half-Maximum and the 
wavelength A: ro « 0.976 p^^^jf^-. There is no simple ana- 
lytical expression for the resulting PSF, and although the 
core bears some resemblance to a Gaussian, the wings are 
significantly largc;r. The atmospheric FWHM varies only 
slowly with wavelength (ro oc A^/^, i.e. FWHM oc A"^/^). 
This results in a negligible dependency of star profiles on 
their spectral characteristics under groimd-based standard 
broadband observing conditions: something essential for 
weak-lensing calibration, as stressed by Kaiser (1999). 

The Aureole 

A so-called "aureole" is observed to dominate all opti- 
cal instrumental PSFs out to distances a few times the 
FWHM away from the centre. This is a feature which 
must be taken into account when simulating deep and 
wide galaxy fields, as it reproduces the background varia- 
tions found on real images around bright stars. Although 
several sources contribute to the presence of the PSF au- 
reole (light scattering caused by aerosols, dust on optics, 
scratches and micro-ripples on optical surfaces, see for in- 
stance Beckers 1995), it is experimentally found to follow 
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Fig. 17. The input pupil used for the simulations. Note the secondary-mirror /primary- focus obstruction and the position angle 
of its support arbitrarily set to 30 deg. with respect to image axes. Left: Pupil function modulus. Right: Real part of the pupil 
function in presence of defocusing, coma and triangular aberrations (PSF_6 from Fig. 13). 



quite a "universal" profile (King, 1971). The intensity 

appears fairly constant too. being close to 16th mag per 
sq.arcsec for a 0th magnitude star. We adopted this value 
here. 



sampling step is that of the PSF (a few mas). Prior to con- 
volution by the PSF, appropriate thresholding is applied 
to the wings of the composite galaxy profile to remove any 
"boxy" limit that may affect shear measurements. 



Pixel footprint 

Finally, the generously sampled (FWHM « 40 sam- 
ples) instrumentaH-atmospheric PSF is convolved with 
the square pixel footprint (0.206", like the CFH12K cam- 
era) and ready to be interpolated over the final image grid. 

Simulated stars 

Stars are simulated within SkyMaker, assuming a constant 
slope of the logarithm of differential number counts of 0.3 
per magnitude interval, and a total sky density of 42000 
per sq.degree down to 7=25. The fairly high slope pro- 
vides a crude match to star counts around 7=20-21 at 
high galactic latitude (e.g. Nonino et al. 1999), while re- 
ducing the density of brighter stars. 



Final image 

A uniform sky background with surface brightness fii = 
19.9 mag.arcsec^^ is added to the image. Surface bright- 
nesses arc then converted to ADUs assuming the image 
is the average of twelve 600s CFH12K exposures, yield- 
ing a magnitude zero-point of 32.6 mag and an equivalent 
conversion factor of 36 e^ / ADU. Poissonian photon white 
noise and Gaussian read-out white noise realisations are 
eventually applied. It is important to note that the simula- 
tion procedure skips three important steps, present in the 
reduction of real images: field warping, image resampling 
and co-addition. These points can have a significant im- 
pact on weak-shear measurements and will be addressed 
in a future paper. 



Simulated galaxies 

SkyMaker galaxies are made of a bulge with an expr"^/'^ 
profile and an exponential disk. All galaxy parameters are 
provided by the Stuff program as described in Sect. 3. 
Obviously, a study of weak-shear systematics requires ac- 
curate profiles, which means that simulated galaxies must 
be well sampled. Except for the most extended ones, the 
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